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INTRODUCTI ON 


This is a semi-annual Status Report for the work accomplished under 
NASA Grant NGR-05-07 1-005 for the Trajectory Analysis and Geodynamics Division, 
Goddard Space Flight Center, National Aeronautics and Space Administration, 

This report contains the mathematical analysis and computation of the 
K=3, order 4; Kt 4, order 6; and K=5, order 7 cyclic methods and the K=5, order 6 
Cowell method and some results of "optimizing" the 3 backpoint cyclic multi- 
step methods for solving ordinary differential equations [2,5, 10] . Cyclic 

methods have the advantage over traditional methods of having higher order for 
a given number of backpoints while at the same time having more free parameters. 
After considering several error sources the primary source for the cyclic 
methods has been isolated. 

The free parameters for three backpoint methods were used to minimize the 
effects of some of these error sources. They now yield more accuracy with the 
same computing time as Cowell’s method on selected problems. 

This work is being extended to the five backpoint methods. The analysis 

I 

and optimization are more difficult here since the matrices are larger and the 
dimension of the '’optimizing space” is larger. Indications are that the primary 
error source can be reduced. This will still leave several parameters free to 
minimize other sources. 
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2. COMPUTATION RESULTS 

Class II Methods on Orbits (Local Correction Error): 

Integrating the two body equations (6.1) for the motion of a satellite 
similar to GEOSB with exact starting values in [ 5 ] we found that with the cyclic 
K=5, order 7 corrector and the PECE algorithm moving from an unstable order 7 
predictor to Stormer order 6 predictor increased the error slightly at small h 
(due to lower order) but greatly decreased the error at larger h (due to in- 
creased stability). This implies that the predictor has a great effect on the 

cyclic method and should be included in future derivations. 

9 

With the latter predictor the PE(CE) algorithm improved the accuracy 
slightly at small h and greatly at larger h. However, iterating to convergence 
with the corrector at each step did not decrease the errors any more than this. 

As can be seen from Figure 1 the cyclic K=5 corrector still needs improve- 
ment as it is being compared to the cyclic K = 4 and Cowell K = 5 with only the 
PECE algorithm. At this h=l sec. the "random" local roundoff error probably 
dominates the local truncation error. The cyclic K=4 error typically oscillates 
in the first few steps while the K~ 5 jumps greatly. Details of Cowell on this 
equation were not available but other results indicate it increases slowly. All 
three soon level off to about the same rate of error growth. 

Class I Methods on y 1 = .5y: 

Because of the strong dependence on predictor and algorithm and because 
the error curves were similar to those for orbits, we decided to being our study 
on linear equations. Figure 2 compares the cyclic (order 9) and Adams (order 9) 
class I, K-5 methods using program CCMPAR in double precision with exact starting 
values. Since 9x4=36 and since the CDC machine carries ^28 places the "random" 
local roundoff error dominates the local truncation error. The h=10 D curves 
are similar to these h=10 ^ curves. 

The Adams coefficients were "perturbed" by 5x10"^^ (see Section 3.1) which 
is the estimated error in the cyclic coefficients. This is probably the reason 
the Adams error greatly increases during the first cycle. 

Class II Methods on y" = y: 

-28 

In order to dominate the truncation and roundoff (^10 °) errors the 

starting eriror =(1,-2 ,4, -8, 16) x 10‘ 20 was propogated using CCMPAR at h=0 in 
double precision for the cyclic, K=4, order 6; cyclic, K-5, order 7; and Cowell, 
K=5, order 6 methods. See Figure 3. 
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Figure 4 shows the results of integrating y" = y with CCMPAR in double 
precision at h=10 - ^ with exact starting values- Since 6x6=36 the roundoff 
error dominates. 

In Figure 5 at h=lO *" with exact starting values the truncation error 

dominates- The cyclic method is better in the first cycle since its order is 

higher but immediately gets worse in the second. 

The cyclic, K=5, order 7-4/5 curves were almost the same as the order 7 

ones so they were not graphed- Single and double precision runs were made 

with h ranging from 0 to 10"^- All methods blew up at 10 ^ (i.e., error after 

- 2 

50 cycles 1) but remained good at h=10 


3. COMPUTATIONAL ANALYSIS 

3.1 Stability With Respect to Perturbations of the Coefficients: 

Because of the propogation of roundoff errors in solving the nonlinear 
stability equations} the K=5 cyclic coefficients were estimated to be in error 
by '—'10 . The K = 6 even more so because of the inaccurate order equation 

solutions. These errors do not effect integrations in single precision 
( ]_4 places) so the conclusions will still hold since they will be based 

in part on single precision integrations. However, the sensitivity of the 
integration errors in double precision to a perturbation of the coefficients 
is a measure of the stability of a method.. 

Table 1 shows the effects of adding to each cyclic, K=5 , order 7 

_ n o 

coefficient 5x10 and to the Cowell, K=5, order 6 ones the same amount 
with exact starting values at h=10~"^ where both roundoff and truncation 
contribute. The Cowell errors increase by a factor of 10 ^ and the cyclic 
by 10~* (less because they were in error to start with). Since the machine 
computes with 28 to 29 places our perturbation is —'10^ times the computation 
accuracy, the computations will now be done 10^ times less accurate, and 
this is about the error increase shown. If this result can be generalized then 
the methods are nearly equally sensitive but not overly sensitive to these 
perturbations- This implies derivations of coefficients must be done more 
accurately than is their practical use in integrating equations- 
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Table 1. Integration Errors 


Error at Cycle # 


#1 


#2 . 

#50 

Cyclic unperturbed 

2 

E-24 

6 

E-23 

2 E-20 

Cyclic perturbed 

4 

£-19 

2 

E- 17 

2 E- 14 

Cowell unperturbed 

2 

E-27 

7 

E-27 

4 E-24 

Cowell perturbed 

2 

£-20 

6 

E-20 

3 E- 17 


3 = 2 Local Truncation Error (Order): 

The order equations were satisfied in all cases and the first few. 
non-zero coefficients in the local truncation error expansions are comparable 
to Cowell's. 


Table 2. Truncation Error Coefficients 





Cycli 

c, K=4 , order 6 



Method # 

= 

1*3 




C 8 

= 

= 0021 

.0026 



c 9 

= 

.0042 

= 0052 



C 10 

= 

.0046 

.0056 



C 11 

= 

.0035 

.0042 





Cyclic. K-5, order 7 



Method # 

= 

1 

2 3 

4 

5 

C 9 

- 

-.0005 

.0026 -.0003 

-.0007 

-.0009 

G 10 

= 

-.0011 

.0064 -=0007 

-=0017 

-.0022 

C 11 

= 

- =0012 

.0085 -.0007 

-.0019 

-.0026 

C 12 

= 

-.0009 

.0079 -.0003 

-.0014 

-.0021 


3.3 Eigenvalues, Vectors, Condition Numbers (Stability): 

The "stability matrix" determines how the local errors are propogated for 

the cyclic methods and for a traditional method used cyclicly. The eigenvalues 

-6 -4 -2 

were of primary concern and were computed- for h-0 and for h-10 ,10 ,10 , 

10“ 1 for the equation y" - y. The h=0 computations were not accurate since 
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the Cowell and cyclic matrices were ill-conditioned, the Jordon block associated 
with A=1 was 2x2, and the computed eigenvalues were overly sensitive to round- 
off errors. The principal condition numbers as derived from the eigenvectors 
(Section 4) were also of importance and were computed at h^lO ' (Section 5.1 
contains more thorough computations in this regard). Finally the row norms of 
the stability matrices varied slowly with h and are: 9.0 for the cyclic K 4 

method, 11.0 for the Cowell K=5 method, and 1723.8 for the cyclic K-5 method. 

The cyclic and Cowell K=5 principal eigenvalues behave the same. The 
Cowell extraneous values remain well within the unit circle at h=10 ^ when it 
blows up in computations implying it is "over stable." The cyclic extraneous 
values leave the unit circle at the same h at which it blows up implying ( 
that the "blow up" point may be moved to larger h by improving the behavior of 
the eigenvalues. 

Table 3. Eigenvalues, A^ 

In the vector (r,i) r is the real part, i the imaginary. 


h 

Root #1 

#2 

#3 

#4 

#5 




Cyclic K=4 



10“ 6 

1.000004 

0.999996 

(-1.0, 1.4E-13) 

(-1.0,-1.4E-13) 

10" 4 

1.00040 

0.99960 

(-1.0, 1.4E-9) 

<-1.0,-1.4E-9) 


io " 2 

1.04 

0.96 

(-1.0, 1 .4E-5) 

(-1.0,-1.4e-5) 


10“ 1 

1.5 

0.67 

(-1.0, 1 .4E-3) 

(-1.0,-1.4E-3) 





Cowell K=5 



0 

1. 

1 . 

-5.E-23 

-7.E-35 

-3.E-35 

10"6 

1.000005 

0.999995 

l.E-28 

-8.E-28 

-l.E-28 

t— * 

o 

i 

4> 

1.00050 

0.99950 

5.E-18 

(-2 . ,-4.)E-18 

(-2. ,4.)E-18 

10“ 2 

1.05 

0.95 

2 .E- 11 

(- 1. , -2 .) E-ll 

(-1. ,2.)E-11 

IO* 1 

1.65 

0.61 

4.E-8 

(-3 . ,-4») E-8 

(-3. ,4,)E-8 




Cyclic K~5 _ 



0 

1 . 

1 . 

-l.E-4 

(0.7, 1.) E-4 

(0.7,-l.)E-4 

10’ 6 

1.000005 

0.999995 

2 . E~6 

(-1. ,2.)E-6 

(-1. »-2.)E-6 

i 

o 

H 

1.00050 

0.99950 

0.00010 

-1.2E-5 

-9.1E-5 

10" 2 

1.05 

0.95 

0.02 

-l.E-5 

-0.005 

10” 1 

1.65 

0.61 

1.56 

-l.E-5 

-0.006 
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Root # 


# 3 


#4 


Cyclic K=4 
Cowell K=5 
Cyclic K-5 


Table 4. Principal Condition Numbers, 


#1 n 

-1.414E-4 +1.414E-4 

6 a 326E-5 -6.324E-5 

-1.033E-8 +1 . 032E-8 


(-0.4,1. 2) 


(-0.4, -1.2) 


4 . MATHEMAT I CAL ANALYSIS AND OPTIMIZATION 


4.1 Summary of Computational Results: 

The Class I and XI, cyclic, K = 5 methods in computations dominated by 
starting error, "random" roundoff error, or truncation error always show 
an error jump in the first or second cycles. The error growth levels off 
to the same rate as Cowell's as you integrate along (increase cycle number) 
on a given equation so this is not a stability problem (Figures 1, 2,3,4 and 
Table 3). The difference between the cyclic and Cowell errors remains a 
factor of about 5xl(P for h=0,10~^ and 10"^ (compare Figures 3,4 and Table 1) . 
At 10“2 (Figure 5) the truncation error dominates the first cycle so the 
cyclic error is smaller but during the second cycle it jumps as before but 

p 

by a smaller factor of about 5x10 . 

In fact it is quite common that the cyclic error actually decreases 
(figures 2 and 3) while this has never been observed for any traditional 
method. This is probably due to the different correctors cancelling the 
errors. It may be possible to choose the correctors so that the errors 
exactly cancel within each cycle thus obtaining an "errorless" method. 

This does not seem to -be possible with any single method. The K = 4 cyclic 
method illustrates something very close to this "ideal" cyclic method. The 
error growth in the first few cycles is smaller and thereafter levels off 
at the same rate as Cowell but at a smaller error level. 

The cyclic method has the largest condition numbers, the smallest 
norm, the slowest growth of extraneous values with h, and the smallest error 
in all commutations. The cyclic K=5 is just the opposite in all respects. 

The Cowell' K=5 is in-between in all respects its condition numbers being 

n t 

about 5xlCT larger than the cyclic at h=10" (Table 4). 
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4. 2 Convergence Proof, Error Bound: 

Some analysis will be provided to explain these results and to 

lead to methods of improving the accuracy. The equation y" = ay 

will be studied in detail and it is expected that improvements will be 

obtained also on orbit problems. 

For y" = f (x , y) the cyclic methods take the form 

(4.2.1) A- Y + A a Y - h 2 (B. Y " + B n Y ") = 0 where 

1 s+1 Os 1 s+1 0 s 

s = 1, 2, ... , S is the cycle number, A^ , are lower triangular 

and Aq , upper triangular K xK matrices consisting of the a^ 

and b. of the K correctors [ 2] , Y g consists of the approximate 

th ' ’ 

solution at the K grid points of the s — cycle, Y are the 

s 

corresponding approximate second derivatives from f (x , y ) , and 

n n 

Y^arethe starting values. 

If each individual method is applied to the exact solution, 
y (x) , restricted to the grid points, we obtain 


(4.2.2) K s 
k=0 


r m m 

U k y(X n+k ) - h b k 


y" (x n+k ) ] = 


t ( x ) 

m n 


where m = 1, 2, ... , K = the number of the method, the local 

truncation error t (x ) = C** 1 y (x ) + C? 1 y l (x ) h + • • • + 

mn 0 7 n I n 


C y^ (x ) h^ + - ■ • [ 9, p. 296 ] . The order of the local truncation 
p n 

error is the power of h in the first non-zero term. The "order" 

of method #m is the order of t -2 . 

m 


Eet T . — C t_ (x __) , t (x ) r 
s+1 1 sK 2 sK+1 


^sK+K-l 1 ] COnSist ° f 
,st 


the truncation errors of the K methods in the s + 1 cycle. Since 
the computer has only finite word length and since we iterate only a 
finite number of times in solving the implicit corrector equation 
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at each step, a roundoff error will be committed. The right side of 

(4.2.2) should be t (x ) + r where r must be considered 

m n m m 


random variables in practice. Let R , = [r_, r_, 

Sri 1 C 


# r ] . 


Writing (4. 2. 2) in matrix form and subtracting (4. 2. 1) we obtain 

(4.2.3) A, E , + A E -h 2 (B_ E" + B n E" ) = T + R . 
v 1 s+1 Os 1 s+1 0 s s+1 s+1 


where E g = [ylx^l - y„ R , 


' * < X sK + K-l> ^sK+K-1 ] a “ d E s 


consists of the exact minus approximate second derivatives. 

For y" = & y (4.2,3) becomes 

(4.2.4) LE s+1 + UE s = T. +1 + R g+1 

2 2 

where L = - q h and U = -cvh or 

(4. 2. 5) E s+1 = A E g + B s+1 where B g+1 = L 1 (T fl+1 + R, +1 ) is the 

-1 - -1 

total local error magnified by L and A = -L U is the "stability 
matrix". In terms of the starting errors, ; 

(4.2.6) E s+1 A E o + s =0 S A B S-s+l . In terms of the 

Jordan cononical form, J, and the similarity transform , P , 

S+ 1 -1 S s -1 

(4.2.7) E s+1 = PJ P E 0+ s = 0 2 PJP B s _ s+1 _ 


(4. 2. 8) Thm: Convergence for y" = ay ; 

If (i) stability: are distinct for h > 0 and | X | 

where 6 ^ ch , 


= 1 + € 


max 


(ii) consistency: the local truncation and roundoff errors || B || £ O (h ) , 

s 

then (iii) the cyclic method converges and the order of the propogated 
error is the min of the orders of the starting error, local roundoff 


error minus 1 or local truncation error minus 1 , and 

3 

(iv) II E s+1 II * K 2 |p. I max [ || E 0 I I + II B. II max (x-x 0 )/Kh]e 


c (x-x)/K 
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r 



where the terms are explained below. 


Proof: Letting be the normalized eigenvector of A associated with 

T T 

\ , Z, that of A , s. = X. Z. are the "condition numbers" , and 
k . k k k k 

p, = 1/s. then the columns of P are X, , I! P I | < K , the rows 
k k k 1 ' 

-1 s -I i s s -1 

of P are p Z [12] , the rows of J P are p k Z , so | | J F | 

K. JK xC xC 

4 

S ^ I P k 4 I max S ^ I P k I max* 1 + 5 )S ' Using inequalitieE for 
the row norm in (4. 2, 7) 

II E II* IIP II II P _1 || [||E 0 || (1 + €) S+1 + ||B || S (1+6)*] 

s = 0 

* I<3/2 K I max t H E 0 II < 1+ch > S+1+ H B s II m ax (S+1) (1+ch)S ] 

■ * k3/2 |P k l m ax (1+ ch) (X ' X ° ,/Kh [ |I E 0 M + (S+Dll B g || max 1 
SK3/2 l P k lmax eC<X ’ X ° )/K[ l |E oH + H B E llmax (x -y /Kh] ' 

S S € 

In the last steps we used (1+ €) <■ e and the fact that if x is in the 

cycle then K(S-l) h £ x - x„ < K Sh . Since | j B |j £ O (h^) the 

second term is O (h) and if the starting errors are at least O (h) then 

we have convergence. End of Proof. 


The stability condition is suggested by Table 3 for a > 0 but 
has not been established yet for a < 0 . It also suggests c = K. This 
convergence proof can be extended to a more general class of f(x, y) by 
using stability at h = 0 (not h > 0) and incorporating h / 0 by using the 
linearization of f (x, y) given by a Lipschitz condition. The theorem 
would resemble that hypothesized in [2] with constants resembling 
those of (4. 2, 8 iv) ♦ Because we are considering linear 
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equations this bound is better than that in [9] for arbitrary f and 
traditional methods in several respects. For example, the order of 
convergence is one higher. Also, for oi > 0 y ~ e X so the 
relative error grows only linearly with x (or with S) as is observed in 
the graphs. 


Both bounds also show if the same local truncation error can be 

obtained with smaller K then the propogated error will be smaller. This is one 

advantage of circumventing the Dahlquist [3] stability criteria in 

addition to that of easier restarting. An advantage of treating the 

cyclic methods in this matrix fashion instead of as an "auxiliary method" 

[4] is that the factor [ p [ is explicit. This is also in agreement 

x max 

with the observations (inc. elliptical orbits) and with Table 4 "explains" 

why methods differ so greatly in the first cycle but all’ level off to 

about the same error growth when the starting and roundoff errors 

dominate. When the truncation dominates the method with the smaller 

truncation error will be best in the first cycle since E - AE + L~\r + T ) 

1 1 10 11 

« L R but in the second E A R R. so the p, factor will enter. 

-l 2 1 k 

This explains Fig. 5. 


Wilkinson [12] presents the theory which shows that the 
eigenvalues of a matrix with larger p will grow faster with respect 

.K. 

to perturbations of the matrix. This explains Table 3 since in our case 

the perturbation is h times the "b" coefficient matrix. Increased 

stability will be obtained with smaller p. 

k 

The matrix approach and particular the error bound with 
p^_ factor explain very well all the graphs and tables containing the 
observed results. 


4. 3 Error in the First Cycle: 

To understand more fully the observed behavior we must study 
the first cycle in more detail than an error bound will allow. At 
h > 0 the are independent so E Q and can be expanded as 
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r 



k = l 2 e k *k and B 1 


Multiplying 


'0 


K _ __ 

k = 1 2 b k X k 


T T 

by Z and using the definition of p gives e = p E n Z , 

K. K ic Kl U K- 

b k = P k B 1 Z k’ and E 1 = k = l 2p k ^k E 0 Z k + B 1 Z k ^ X k 


A similar expression could be obtained for E from (4, 2. 6) 

S ^ 

in which would appear implying that only the principal condition 

numbers matter . It may be possible to pick such that the 

expression in brackets is 0 . In the general case there are two 

difficulties with this approach: (i) if truncation error £ roundoff 

error the second term is random and there is no hope of estimating 

it, (ii) if truncation > > roundoff error the second term will be 

extremely difficult to estimate requiring knowledge of the higher 

derivatives of the solution. If this is the case it may be possible to 

go to a more accurate method at the same h (higher order, smaller 

C. / 0 , or a completely different method) . The '‘cancellation" 

effect (4. 1) in the cyclic methods could be taken advantage of in either 

(i) or (ii) . 

K ST 

Supposing the local errors are < < E^;E^ = ^ \ ^^k^k ^O^k 

This also represents the propogated effect of the local error at a single 

step. For simplicity we will study E^ = (0, 0, 0, 0,1) for the Cowell 

-4 

and cyclic K = 5 methods at h = 10 using Tables 1, 3, and 4. The 
products of the extraneous p X were < . 01 of the principle product 
even in the first cycle so they will be ignored. The first componant of 
E for Cowell = (6. 326 x 10' 5 ) _1 (1. 00050) (. 70714) (. 4471) 

+ (-6. 324 x 10’ 5 )- 1 (0. 99950) (. 70707) (. 4473) 

= (1 + 5h) (. 4471) (. 70707 + . 7h)/(. 6324h) (1 + -|t h) 

t o«5 

- (1 - 5h) (. 4471 + 2h) (. 70707)/ (. 6324h) 

» (.4471) [.70707 + (5 x . 7 + . 7) h + 3. 5h 2 ] [1 - h ] /. 6324h 

• t) J 

- (. 70707) [.4471 + { 2 - 5 x . 44) h - 10 h 2 ] / .6324 h 



« [( . 4>(. 7)(- -^-h) + .4(5x . 7 + . 7) h (1- -|rh) - . 7(2-5x. 44)h]/. 63h 

• Dj ■ Qj 


& 7. We see that although s.. ~0(h) the fact that 

\\ I- 1^1 ~ l x ! I - |X 2 |~|Z 1 | - \Z 2 | ~ O(h) , that 

2 

| s | - | s | ™ O (h ) , and that there was an overall difference of sign 

X Lr 

all lead to the cancellation of the zeroth order term in the numerator 
thus making the error reasonable; a factor of 7 is in agreement with 
Table 1. 


It may be possible to apply perturbation theory [12] to verify 
the above conditions (perhaps even the sign difference) on arbitrary 
f(x, y) , provided the p^_ are not too large# since X^ = X^ and 

X = X at h = 0 , The X and Z pairs must have nearly - 

12 

corresponding componants, 

For the cyclic method | X^ | - | X^ | ^ | X | - I ^2 I 


Z 2 ~ 


z 

O (h) # s, | - s ~ O (h ) , and there is a net 


difference in sign so also here 

(4. 3.1) | | | | ^ c h || Eq || / | s^ | where c = constant + O (h) 

includes a term | + s^ | / | s^ | h so E^ includes a term 

2 

l|E n || |«, + s | / | s | * For the cyclic however 

v 1 C$ J, * 

J | E^ | j ^ c 10 4 /10 ^ 1 1 Eq 1 1 || E^ | j x 10 4 which is in surprising 

agreement with Table 1 in view of the simplifying assumptions made. 


4.4 Optimization Criteria and Methods: 

The problem facing us in not one of using some parameters to 
satisfy a set of linear or nonlinear equations as almost all procedures 
for deriving numerical methods are. We do not know what values can be 



attained by criteria such as jj A j| nor do we know how many 
parameters it would take to solve such an equation. The former 
problem is surmountable by just assuming smaller and smaller 
values from the present one. However, the wrong choice of number 
or type of parameters would make an equation solution impossible. 

^ Although nonlinear optimization procedures are computationally 
more complex and lengthy than solving nonlinear equations they do 
solve the above problems. There are many ways to state our problem: 

(i) minimize il A|| subject to order and extraneous eigenvalue 
conditions, (ii) minimize extraneous eigenvalues subject to order and 
l|A|| conditions , (iii) minimize the | p | factors subject to order 
and eigenvalue conditions or vica versa. The condition of A or the 
norm of JL or other ways of stating the primary error source could 
be included. Auxiliary conditions could include local truncation error, 
stability of the eigenvalues, the size of the coefficients themselves, 
and perhaps even some conditions on the behavior with respect to 
random local error (roundoff) . For K = 3 a simple mapping program 
was used to ’'optimize 1 ' certain criteria. This is too expensive a process 
if reasonable accuracy is required at K = 3 and for any accuracy at K = 5 . 

The iterative, nonlinear optimization algorithms seem to fall 
into three classes: (i) require no derivatives , (ii) require the 

gradient of the object function, and (iii) require second partials of the 
object function. Methods of (i) are slow converging and require as 
many function evaluations as (ii) [11] . In our problem we do not have 
an explicit expression for the object function much less its derivatives. 
Difference quotients have been used in (ii) but would be very inaccurate 
in (iii) [7] . The programming of methods (iii) is also complex. The (ii) 
methods, then, consist of the linearly convergent steepest descent (in 
the negative gradient direction) [1] and the quadratically convergent 
conjugate gradient [8] methods. 
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Several modifications of the latter are being tried under the program 

name OPTIMA. We are trying to reduce the total number of function 

evaluations since these will probably involve computation of eigenvalues 

and vectors which is very expensive. For example, to optimize a 

K = 5 method using 5 parameters will require ^ 12 eigenvector 

computations per iteration for, perhaps, 200 iterations « 2500 

evaluations at about 1 second each at about $.20 each second ^ $500; 

* 

which does not satisfy budget constraints* 

Once this program is working the constraints must be added. 

The best way to do this seems to be the penalty function method. 
Whenever the minimum search wanders outside the region where 
constraints are satisfied a penalty proportional to the size of the 
constraint is added to the function we are trying to minimize [6] . 

This tends to keep the search within the constraint boundaries. 

The procedure will then be to solve the order equations 
parametrically, minimize | | A 1 1 subject to extraneous eigenvalues 
^ say *“ . We will then work on the more expensive minimize j p . | 

procedure which promises greater improvement. If imporvement up to 
Cowell is obtained then other optimization criteria can be added. If 
no improvement, order will be dropped to Cowell's and the above 
repeated. At this point, this work would begin to blend with optimization 
of traditional methods [5] so if no improvement is obtained over 
the already optimized traditional methods, then this phase of the work 
will be terminated. 


5. OPTIMIZATION RESULTS 

Using OPTK3 several K = 3 order 4 methods were derived 
that are both better and worse than Cowell with respect to | ] A 1 1 , 

| s^ | and I s^ + s 2 I afc various h on y" = y (Table 5), The methods 
are not optimal due to the crudeness of the program, however improved 
computational accuracy is shown for some of the methods. 
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|| A || did not vary much over the h interval considered. Also 


S 1 = 1/d l 

, the larger 

the better, 

and | s, + s„ 
1 1 2 

} should be ^ 

• o(h z ). 

the smaller the better 

, The integrations were 

done in double precision. 

The errors are 

the last componant at the first and 100th cycle 

S. 

Method #7 

is Cowell's 

* 






Table 5 

Some K 

= 3, Order 4 Methods 


Method 

llA|| 

h 

1 s x 1 

1 s i + s 2 1 

error 1 

error 100 

# 







1 

9.0 

0 

4. 0E-15 

5. 6E-26 

5E-20 

4E-17 



E-6 

2. 7E-7 

1. 7E-12 

IE - 27 

IE-24 



E-4 

2. 6E-5 

2. E-8 

IE- 25 

6.1E-22 

2 

9. 0 

0 * 

1.4E-14- 

9. 2E-28 

2E-19 

2E-18 



E-6 

1. 2E-6 

1. IE -11 

6E-28 

3E-25 



E-4 



8E-26 

6. 2E-22 

3 ' 

4. 8 

0 

3. E-14 

3. E- 27 

6E-20 

IE-17 



E-6 

6. 8E-7 

4. 3E-12 

2E-27 

2E-23 



E-4 

6. 8E-5 

4. 3E- 8 . 

IE-25 

6. IE-22 

4 

6. 2 

0 

1. E-14 

1. E- 28 

2E-19 

3E-18 



E-6 

1. IE- 6 

8. IE -12 

2E-28 

2E-25 



E-4 

1. IE-4 

8. IE- 8 

9E-26 

6. 3E-22 

5 

6. 7 

0 

8. E-15 

2. E-27 

2E-19 

2E-17 



E-6 

1. 0E-6 

5. 6E-14 

2E-28 

2E-24 



E-4 

1. OE-5 

5. 6E-10 

9E-26 

6. 5E-22 

6 

6. 4 

0 

1. E-14 

2. E-27 

2E-19 

IE -17 



E-6 

1. IE-6 

6. 5E-13 

7E-28 

6E-24 



E-4 

1. IE-4 

6. 5E-9 

9E-26 

6. 4E-22 

7 

7. 0 

0 

9. 9E-15 

1. oE-28 

2E-19 

2E-17 



E-6 

8. 2E- 7 

8. 2E-13 

3E-27 

IE-23 

Cowell 


E-4 

8. 2E-5 

8. 2E-9 

8E-26 

6. 2E-22 
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The spread between values is greatest at h = 10 


- 6 


as is 


the spread in computation errors. Methods with larger s are better. 

- 6 ^ 

For example methods #2, 4 at h = 10 have the largest s^. and the 
smallest errors being a factor of 50 better than Cowell at cycle 100. 
Method #1 is an exception having worse s and s + s than Cowell 

^ J. X M 

but smaller error at h = 10 . Methods that are ‘'better" at one h 

seem to maintain this at other h also. The computation errors at 
- 2 

h = 10 were all the same perhaps due to dominating truncation errors 
or lessening effect of larger s^ . 

Eigen computations with methods #1, 4, 6, and 7 show 

extraneous values remain at 0 and the principal ones = 1 + 3 h as 

( 

expected. The componants of the vectors X^, X^, Z^, and Z^ change 

as O (h) with h so the analysis of (4, 3. 1) will apply. 

These results imply improvement in the condition numbers 
can give smaller integration errors. 


6. SOFTWARE DESCRIPTIONS 

ELLIPSE: 

‘ Integrates two body orbits with Class II cyclic or traditional 

predictor-corrector methods. The major part of the program is the 
same as the GEOSTAR subroutine CSTEP documented in [2 ] . The 
program uses exact starting values and prints the error in each of the 
three position components by using fixed point iteration to solve 
Kepler's equation. It has been modified to correct more than once. 
Computations were in double precision (about 16 places) on the IBM 360 . 

COMPAR: 

Integrates y ! = a y with Class I cyclic or traditional 
correctors and y" = ay for Class II. The starting values are taken 
backwards from x = 1, 2 so that methods of all K start integrating at 
the same x value allowing direct comparisons in the first cycle. 
Starting values and the error at each step are found using the analytic 
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solutions. Given the coefficients a^ and b^* 1 the stability matrix, 

A, is formed as in (4. 2. 5) depending on the Glass and a . It is applied 

once each cycle to obtain the solution. It is possible to add a 

specified local error, B , in each cycle in addition to the normal 

s 

roundoff and truncation errors. If h = 0 the process implied by 
(4.2.6) is undergone where E^ and can be specified. They 

were usually chosen to be (1, -2, 4, -8, 16) E-20 and 0 respectively, 

however, at times B = E n . Computations were in single (about 14 

s u 

places) or double (about 28 places) precision on the KRONOS time 
sharing system on the CDC 6400 . 

EIGENP: 

Computes eigenvalues and vectors for each stability matrix, 

A, using the QR, double step, inverse iteration algorithm [12] . The 

Class II matrices at h = 0 were ill-conditioned so the program either 

did not work at all or gave poor accuracy. Program HESSEN (QRIEG) 

supplied by Mel Velez was used to verify (eigenvalues only) some cases. 

A later version of the program, EIGEN, also computes eigenvectors 
T 

for A and the condition numbers . Computations were done in extended 
double precision (about 32 places) on the IBM 370 and in double precision 
(28) on the KRONOS system. An input parameter specifies the 
approximate number of places of desired accuracy before iteration 
terminates. We usually got by with 20 decimal places. Computing 
time for a 5 x 5 run with EIGEN « 1 sec. & $, 20 on KRONOS. 

OPTK 3 

Computes II A || , max extraneous root, | | , and 

| s^ + | for a specified h for y u = ay using the parametric order 

equation solutions at each grid point of an input specified grid for the 
free parameters . The program functions for Class I or II at any K or 


- 18 - 



order when the correct order equation solutions are given. Using the 
program in a man-machine, time sharing mode the user specifies 
the .initial grid by typing in the left and right hand endpoints and step- 
size for each of the free parameters. After visually inspecting the 
output for the "best" region he immediately specifies a finer grid in 
this region and repeats, obtaining better methods. EIGEN is called at 
each grid point so this program was too expensive at K = 5 costing about 
$10 per run for very coarse grids. The KRONOS system was used. 

OPTIMA: 

Still being developed. Minimizes a nonlinear function of n 

variables, f, for which no explicit form is given. In the ith iteration 

one must choose a stepsize to approximate the partial derivatives of 

f, A x , a downhill direction, D , and a step length in this direction, 

z. After considering several algorithms [1, 6, 7, 8, 11] it was 

decided to start with the method of steepest descent where D. = - Vf(X.) 

and z was found by minimizing a parabola fit thru X., the directional 

i 

derivative in the direction ( = slope of parabola), and a second, 
arbitrarily chosen point in the D. direction. Per iteration, this 
algorithm required only n + 2 f evaluations. Convergence was good 
on simple quadratic f for small n but slow for the function with a narrow, 
curved, ’’banana' 1 shaped valley [8], For f = 1 1 A | | for K = 5 , n = 10 
with no constraints the parabola minimum overshot f min. so often 
that Xj, was changing too radically to converge. The problem here 
is that || A || is a high degree polynomial in the free parameters, too 
steep for a quadratic. Even if we could get down these steep walls we 
would be in a narrow, curving valley like the "banana". 

The program was modified so that D. = - Vf (XJ + 

Eh ^ | V f PC ) | / | V f (X. | ^ called the conjugate gradient 

direction [8] . Theory states that if f min lies in a long, narrow. 
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quadratic valley then convergence is assured in ^ n iterations while 
the steepest descent will take many more. For general functions one 
must set D. = - V f every n + 1 iterations. This algorithm did converge 
faster on the banana but still overshot on j | A jj . Instead of a parabola 
fit, Fletcher uses a cubic fit to two points in the D. direction and their 
directional derivatives. This requires ^ 2nt 2 f evaluations and 
possibly 2n+2 more since he sometimes fits another cubic to obtain 
a better approximation to f min in the D direction. Our algorithm is 
being modified to successively fit parabolas to obtain a better 
approximation to f min in the D. direction. We have taken A x = \ z 
and this seems to work. 

With the latest version of our algorithm after 5 0 iterations 
with 50 (2 x2) = 200 f evaluations we converge to the same point 
on the banana as Fletcher does after 20 iterations with 
s 20 (2x2 + 2) ^ 120 f evaluations. TIis actual number of f evaluations 
could be as high as 200 but he does not give this data. 

Part of the programming support for ELLIPSE and OPTIMA was 
cosponsored by the NASA-PSS contract, part of EIGENP by C, Shipp 
for graduate course credit at CSUF, and the remaining by E, Spiehler 
under the NASA- CSUF Grant, 
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Figure 4. : Cyclic vs. Cowell on y n = y at h=10’^. Roundoff error dominates. 
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